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In this article we present a caracterization of the vibrational spectrum of the H^" molecule using 
the correlation function quantum Monte Carlo (CFQMC) method and a genetic algorithm study of 
the topology of the potential energy surface used in this work. The vibrational modes associated 
with the HjJ"-H2 torsion and stretching posses very flat minima. As a consequence the fundamental 
frequencies corresponding to these modes are poorly described in the harmonic approximation. The 
vibrational frequencies obtained in this work are in good agreement with the available experimental 
data as well as other computational methods found in literature. In our genetic algorithm study of 
the potential energy surface using cartesian coordinates we have found some unexpected minima. A 
careful analysis shows that some of these minima are described by the same curviliniar coordinates 
in which the potential is described. However, they represent nonequivalent molecular geometries. 

I. INTRODUCTION 

The presence of a molecular or atomic ion in a H2 atmosphere results in the formation of hydrogen ionic clusters 
X + (H 2 )„. Properties of these clusters, such as solvation-shell distributions 0, Q, binding energies 0, 0, [f|, and 
vibrational spectrum 0, |f| , have been theoretically investigated for a large variety of ionic cores X + . 

The H2 units are connected to the core mainly by means of monopole/induced-dipole interactions, although there is 
some charge transfer from the H2 units in the first solvation shell to thepositive center. As a consequence of this kind 
of interaction, binding energies lie on the range of 3-10 kcal/mol 01113 that is not enough to destroy the molecular 
identity of each H2 unit. On the other hand, these binding energies are enough large to make the hydrogen clusters 
potentially useful in storage of hydro gen fuel, by taking advantage of the cluster formation in the positive-charged 
sites of adequate crystalline materials |lCt ITU . 

The first species of the homogeneous series H^~(H2)„, the H^j", is an atypical hydrogen cluster ^2- Its electronic 
properties does not fit in those of the large members of the series 0, @ . As Ohta et al. have recently shown , one 
main reason for such anomalous behavior is the isomerization process 

H+(H 2 ) ^ H 2 (H+)H 2 ^ (H 2 )H+, (1) 

which take places even for very low temperatures. 

Since its detection in 1962 jl3| . many experimental studies of dissociation energies and thermochemical properties of 
the H+ clusters have been performed 0, Hj, UM ■ For the purpose of the present work, the investigations performed 
by Okumura et al |14| using infrared vibrational predissociation spectroscopy are particularly relevant. They observed 
the broad and structureless vibrational bands corresponding to the H2 and the H^" units into the (H. 2 ) n ionic clusters 
for n = 1-6 and analyzed the dependence of the shift of these vibrations with the size of the cluster. A similar approach 
was latter used by Bae [l7j . 

Numerous ab initio calculations have been carried out 0, 0, S 0, IM US, E3, EM EM E3, H3, H3, |2^| 

aiming to determine their equilibrium structures, low-lying stationary points and dissociation energies. A very detailed 
ab initio investigation of the H^" cluster was carried out by Yamaguchi et al !3jj. They revealed a quite complicated 
structure of the interaction potential: ten distinct stationary points were located. Their high-level calculations 
predicted a C 2v structure as the global minimum. Spirko et al. studied the potential energy surface of and 
its corresponding vibrational spectrum using the multi-reference configuration interaction (MR-CI) methodology 
|25l 1261 127| . Besides these works there are only a few theoretical papers dedicated to the infrared spectrum of and 
they are all retricted to the harmonic frequencies 0, H, ITTl l2ll |22| . A recent work by Barbatti and Chaer Nascimento 
examined in some depth the vibrational spectra of H^ n+1 using the VSCF methodology |fj. As mentioned above, in 
addition to the scarce theoretical data only the experimental results of Okumura EH an d Bae ^3 are available for 
the vibrational spectrum of this system. Therefore a more complete caracterization of this system is warranted. The 
focus of this work is to provide a detailed description of the vibrational spectrum of . 

A full description of the vibrational spectrum of H^" needs to take into account nine degrees of freedom. The 
computational effort in the traditional variational methods that rely on integration using normal quadratures is 
prohibitive. One of the alternatives to overcome such shortcoming is the use of Monte Carlo methods to calculate the 
multidimensional integrals which appear in time independent problems, |28l I29I l30l l3ll I32I l33| which do not present 
the same limitations of the variational methods. This method has been applied successfully for the computation of 
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the vibrational and ro- vibrational spectra of triatomic and tetratomic molecules [3(], EI HH l3~H l35l l36l l37l l38| . In 
this work we used the CFQMC methodology to obtain the vibrational states of the H^" molecule. We also study the 
PES used in this work using the genetic algorithm|ii|H3,|EllEllEllEl. This article is organized as follows. A short 
description of CFQMC methodology and of the genetic algorithm is given in the next section. Sec. Illll is devoted to 
the presentation of the results, and in sec. IIVI the discussion and concluding remarks are presented. 



II. METHODOLOGY 



A. Correlation Function Quantum Monte Carlo 

We start from the Born-Oppenheimer nuclear Hamiltonian of a molecular system (in the center-of-mass reference 
frame) 

# = -E^- v ' + ^ R )' ( 2 ) 

i=i ^ l 

where R is the vector of the coordinates of all particles of the system, 1^(R) is the potential energy of the nuclei. The 
vibrational energy levels of the system will be then computed by solving the eigenvalue problem: 

H&i(B.) = £7j$j(R), (3) 

where Ei and $j(R) are the eigenvalues and eigenvectors of H. Given a trial basis set of m known functions /i(R) 
we can define the following generalized eigenvalue problem 

m 

J^lHij - AkN^dkj =0, (4) 
i=i 

where dk is the k th eigenvector and its associated eigenvalue, and 

N l3 = J dR/i(R)/ :) (R), 
Hij = / dR/ i (R) J E/j(R) , (5) 



are the overlap and Hamiltonian matrix elements associated with the basis set. The matrix elements defined in 
Eq.© are evaluated using Monte Carlo integration techniques. The only difference of this method and the traditional 
variational methods is the use of Monte Carlo to compute the integrals. Of course an accurate spectra will depend on 
the quality of the basis set {/^(R)}. In this work we use the basis set of reference [3(3]. Namely, for the ground state 

il>o = ex P C£2 ^S^A^AS^) (6) 

with AS„ = S„ — 5°, and = — rj\ the distance between atoms i and j, and S® — — rj\ the equilibrium 
distance between atoms i and j. The variational parameters A ufl are optimized in order to minimize the variational 
energy or its variance. The trial functions for the excited states are given by 



This basis set was app lied successfully in the study of vibrational spectra of 2-, 3- and 4-atom molecules 28, 30, 31, 
l32l I33L I34I I35L l36l l37l l3S| . An improvement to this methodology is the use of the diffusion Monte Carlo techniques 
to "project" out the excited state spectrum. We restricted ourselves to the variational implementation of the method 
in order to have a comparison in the same footing of this technique with the VSCF, which is also a variational based 
method. 

Another important aspect related to the accuracy of the CFQMC is the quality of the potential energy surface 
(PES) describing the motion of the nuclei. In this work the PES of Spirko et al. |2g fitted to MR-CI all-electron 
computations at 110 different configurations of the nuclei was used. A set of curvilinear coordinates was used to 
facilitate the computation of the vibrational spectrum . The purpose of the CFQMC calculations is twofold. First, 
to show that the method can be applied to compute the spectrum of a 5-atom molecule. Second, to determine the 
accuracy of the PES of ref . j2(| . 
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B. Genetic Algorithm 

In our genetic algorithm the population is coded in a binary discrete cube named A, with I x m x n bits. The 
elements of A, a^, are either or 1, with i,j, k integers numbers 1 < i < I, I < j < m, 1 < k < n. The label i refers 
to the component i of the gene j of the individual k. Therefore, A represents a population of n individuals, each one 
of them have a genetic code with m genes. Each gene is a binary string with I bits. 

The genetic code of the individual k is given by 

[a] fe = [aik,azk, —,a m k], 

were 

l 

= X] 2 t_1 Oyfc (8) 
i=i 

is a integer number composed with the binary string aijkci2jk--aijk--aijk- It is defined on the interval [0,2 l — 1]. To 
define the real search space for each parameter, we transform 

a jk -> a k = ^— a jk + cij (9) 

where ajk is a real number defined on the interval Sj = [a™ m , a™ ax \. 

We define the phenotype of the individual k, Vk = V'QaJk) where [a]k = [aik, &2kj a jk, a mk] is a set of 
coefficients that characterize the individual k. With this we define the fitness of a phenotype k (set of coordinates) 

Fk — S max Sk 

where Sk is the energy given by the SEP for this phenotype, and S max is the worst individual in the population. 

We use the most common operators: selection, recombination and mutation. The selection operator normalize the 
vector Sk 

p >=ik < 10 » 

that represents the probability of each individual been selected for a recombination through a roulette spinning. For 
the purpose of this work we selected n/2 individuals (parents) that will generate, through the recombination operator, 
n/2 new individuals (offsprings). So, to make a new generation we joint the n/2 old strings (parents) with n/2 new 
strings (offsprings) in order to maintain the population with fixed number n. The recombination operator is a cross- 
over operator that recombine the binary string of each gene j of two random selected individuals to form two new 
individuals. In this work we use a two random point cross-over. 

The mutation operator flip N mut random selected bits in a population. We choose N mut to make the probability 
of change of a given bit equal to 0.01 per cent. So, in a population of I x m x n bits, we make 



9=1 11 

I x m x n 



where q is the probability of change of one bit. 

An elitist strategy is used. It consists of copying an arbitrary number N e i of the best individual on the population 
in the next generation. It warrants that this individual will not be extinguished. 

We found a large number of acceptable solutions. The set of all solutions is the definition of search space (r). 
The length of T is defined by the number m of coefficients and the length I of the binary codification. Each one of 
the 12 coefficients (coordinates), that define the individual k, can assume 2 l distinct values. So, an individual in the 
population is only one possibility among 2 lxm . This value defines the length of T. The length of T describes the 
number of digits that are used to express a real value ajk and shows the minimal difference between two possible 
values of ajk- Being each coefficient defined on an arbitrary interval Sj, the precision of the coefficient ajk is 



r max r 7nin 

U j U 3 



2' 



(12) 
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If we do not have any information about the order of magnitude of the ajk values, we must choose the 6/s intervals 
such that they cover the greatest number of values. However, after some generations, we obtained more precise 
information about the order of magnitude of each coordinate ajk- In order to improve the performance of a standard 
genetic algorithm (GA), we include in our technique the concept of dynamic search space. It consist in the use of 
information of past generations to determine the length and precision of the search space for the next generations. For 
the first generations, when we have few information about the magnitude of the coordinates, we do not need many 
digits to represent a real number ajk, that is, we use a low precision codification given by a low value of I. In this way, 
we make T a "small" search space and the GA can find the regions of acceptable solutions faster. Once found some of 
these regions we can redefine the Sj's intervals and rise the precision rising the length of binary codification I. After 
extensive trials of the parameters values we take m = 12, n = 100, q = 0, 01 and N e i — 10. Beside that, we always 
start the GA with a random population defined in the initial intervals Sj — [a™ 71 , a™ 11 *] = [—10, 10] and set the initial 
value for the length of the binary codification I = 12. In this way we had a search space of length 2 lxm = 2 12x12 = 2 144 
and the minimal difference of two possible values of ajk is ^rr — 49 x 10~ 4 . After 1000 generations we redefine I = I + 4 
and 5j = [a™ in , a™ ax ] where a™" = a ]best +a™ m x 10" 1 , a™ ax = a ]hest +a™ ax x 10" 1 and a jbest is the fitest individual 
in the population found along the last 1000 generations. We set 10000 generations for each run of the GA. It should 
be pointed out that the algorithm is very robust and works properly with an wide range of these parameters. 

III. RESULTS 

For well-behaved systems a first, usually good, description of the fundamental vibrational frequencies is the harmonic 
approximation. In Table I we present the harmonic frequencies corresponding to the PES of ref. pr| . used in our 
CFQMC computations, together with those obtained at MP2/6-311G(d,3p) level as the first step in the VSCF 
methodology, the ones obtained by Prosmiti et al. using the QCISD(T)/cc-PVQZ [^j, and the CCSD(T)/aug- 
cc-pVTZ of Prosmiti et al 5J. The results are in reasonable agreement with each other. The largest discrepancy 
is observed in the normal mode corresponding to the second skeletal motion. Whereas the value obtained using 
the PES of ref. [U is 1660 cm" 1 , the frequencies obtained using the MP2-VSCF, the QCISD(T)/cc-PVQZ and 
the CCSD(T)/aug-cc-pVTZ are 1201, 1170, and 1174 cm -1 , respectively. This difference may be due to either the 
difference in the levels of computation, or to the fitting of the PES itself. Unfortunately, in ref. j2(J the frequencies 
in the harmonic approximation are not presented and this question may not be uniquely answered. To gauge the 
validity of the harmonic approximation one must include the anharmonic effects. These effects were considered at the 
CFQMC level of computation. 

In Table II we consider the anharmonic effects at the CFQMC(VMC) level of calculation. For comparison, the 
results of Spirko et al. [26| methodology and some available experimental results are also included. As mentioned 
above, the computation of the spectra was performed using the CFQMC method in the variational form. The basis 
set consisted of 220 basis functions and the parameters of the wave function used in our calculations were adjusted 
to minimize the energy of the ground state and the first nine excited states. As one can see, there is a reasonable 
agreement between all the computations. In particular, the results for the 8 th and 9 th states obtained in reference 
|26|. which correspond to the Hg~ symmetric stretching and the H2 stretching motions, are in good agreement with the 
experimental results of ref. ^4| . The CFQMC results differ a little and they are closer to the values of the frequencies 
of the H2 and H^~ isolated molecules. Although the CFQMC of these frequencies in the H^" are shifted in the right 
direction the shifts are much smaller than the ones predicted by the experiment. 

In FIC^wc display the fundamental frequencies of the H2, H^", and H^" molecules in the harmonic approximation, 
using the CFQMC, and the experimental results [l4|- A similar analysis of the MP2-VSCF results is presented 
in ref. In all cases the frequency corresponding to the symmetric stretching of H^" is blue-shifted in the H^" 

complex and the one corresponding to the H2 stretching motion is red-shifted. However, the shifts predicted by 
the CFQMC are much smaller than the ones predicted by the harmonic approximation, the experiment and the 
MP2-VSCF computations. This difference may be due to the coordinates used to fit the PES used in the CFQMC 
computations. This PES uses a set of curvilinear coordinates |26j which are ill-defined as two completely different 
and non-equivalent set of configurations of the five atoms correspond to the same set of curvilinear coordinates. In 
the case of the harmonic approximation or in the original calculations of ref. pr| , the computations are performed 
close to the minimum and this ill-definition of the coordinates do not affect the final frequencies. In the case of the 
CFQMC the integral is evaluated in the whole space, therefore, making the results more susceptible to the choice of 
coordinates used to describe the PES. Nevertheless, our calculations seem to be describing the stretching motion of 
the H2 and the symmetric stretching of H^" as independent of all the other vibrational modes. 

The CFQMC results are limited by the quality of the PES but has the advantage that one can compute the 
full spectrum, limited of course by the size of the basis set. In FIG. |3we present the results of all the vibrational 
frequencies obtained in our CFQMC up to 5000 cm -1 . The lines are broadened by 30cm _1 . The full arrows indicate the 
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fundamental frequencies as obtained by our CFQMC computations, while the dotted arrows indicate the experimental 
results |14|. The experimental results fall in overtones of our computed results rather than fundamental frequencies. 

Because of the disagreement between the CFQMC computed frequencies of the H2 stretch and the symmetric 
stretching H^" in H^" , with the experiment and other computations we decided to perform an additional study of the 
topology of the PES used in this work. We decided to search for other minima besides the accepted C2V minimum. 
This search is performed with the genetic algorithm which has been demonstrated to be a very robust method for 
global minimum search. 

The inital population was completely random and after a few generations we observed the different minima (local 
and global). Fig. [3^ shows two candidates for the global minimum. The energies of the two minima are identical. 
Further analysis shows that in the curvilinear coordinates in which the potential is defined they are identical minima. 

Based on the minima of Fig. [3^ we started a new population and run the program for a few more generations. The 
minima of Fig. are again among the lowest energy structures. However, a new lower energy minimum is found. 
This minimum is displayed on Fig. As one can see this candidate for global minimum is a configuration in which 
the H2 and the H^" molecules appear as not bound. This explains why the CFQMC computations frequencies of 
the H2 stretch and the symmetric stretching Hjj~ modes are very close to their free values. The CFMQC integration 
samples configuration in the whole space and they will include the structure of Fig. \3]p. The same does not happen 
with the harmonic frequencies as they were computed near the accepted C2V global minimum. 

IV. CONCLUSION 

We have obtained the vibrational energies of the cluster in the harmonic approximation and by using the 
correlation function quantum Monte Carlo (CFQMC) methodology. The lowest frequencies of the H^" are not well 
described in the harmonic approximation due to the flatness of the potential energy surface of the complex. The only 
available experimental data for this system are the frequencies of the stretching motion of the H2 and the symmetric 
stretch of obtained by IR predissociation experiments in the H^" conplex. These frequencies are red- and blue- 
shifted when compared to these frequencies in the free molecules. We have seen that the shifts using CFQMC are 
smaller than the shifts obtained by the harmonic approximation and the experiment. 

Our genetic algorithm search of a global minimum revealed two interesting facts about the PES of ref. [2(j. First, 
their choice of curvilinear coordinates is ill-defined, non-equivalent configuration in cartesian coordinate have identical 
curvilinear coordinates. Second, and most important, the global minimum is a configuration in which the H2 and 
molecules are not bound. This explains the small values of the shifts of the frequencies of the stretching motion of 
the H2 and the symmetric stretch of H^" when compared to their values in the free molecules. 

In conclusion, our CFQMC and genetic algorith study of the PES of ref. |2(j shows a good agreement with other 
computations and the available experimental data. The differences are explained in terms of the definition of the 
curvilinear coordinates used in the definition of the PES. As a future work we are proposing to reparametrize the 
PES in more appropriate coordinates with additional CFQMC computations. 
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FIG. 1: Comparison between the frequencies of the H2 and Hjj~ molecules within rlj cluster with H2 and H^" isolated molecules, 
a) Harmonic approximation to PES of ref . |26j| . b) CFQMC computations using the PES of ref.|26j|. c) Experiment |T^|. Full 
lines represent the frequencies of the H^J" cluster. Dotted lines represent the frequency of the isolated H2 molecule. Dashed 
lines represent the frequencies of the isolated molecule. 
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FIG. 2: The full CFQMC-computed vibrational spectrum of the cluster up to 5000 cm 1 . The lines were broadened by 30 
cm -1 . The full arrows correspond to the fundamental frequencies. The dotted arrows correspond to the measured frequencies 
of the H 2 and in the H;j" cluster 

TABLE I: Vibrational frequencies in the harmonic approximation of the Hg" (cm -1 ) 



Mode 


PES of Ref. 


[26] a MP2" 


Prosmiti I c 


Prosmiti II d 


Main feature 


1 (aa) 


215 


211 


206 


206 


Hj-Ha torsion 


2 (ai) 


477 


482 


495 


502 


H^-Iib stretch 


3 (6a) 


848 


822 


812 


815 


Hjj" rocking 


4 (fax) 


970 


881 


866 


868 


Skeletal motion 


5 (60 


1660 


1201 


1170 


1174 


Skeletal motion 


6 (ai) 


1768 


1862 


1838 


1840 


IlJ bending 


7 (fe.) 


1819 


2186 


2131 


2134 


H^" asym. stretch 


8 (ai) 


3714 


3761 


3668 


3670 


sym. stretch 


9 (oi) 


4177 


4247 


4115 


4118 


H2 stretching 



a Present work. 



b MP2/6-311G(p,3d), present work. 
c CCSD(T)/aug-cc-pVTZ [gl. 
d QCISD(T)/cc-pVQZ Ref. Q- 
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FIG. 3: a) Two candidates for global minima of Hj" molecules, b) Global minimum of the H^" as denned by the PES of ref. 

m 



TABLE II: Vibrational frequencies of the H^(cm ) 



States CFQMC" 


Spirko f26 


MP2-VSCF a 


Exp [14] 


1 598 


622 


596 




2 676 




776 




3 817 




879 




4 952 


973 


989 




5 1058 


1238 


1142 




6 1471 


1383 


1173 




7 1738 




1735 




8 3281 


3471 


3515 


3532 


9 4117 


3897 


3921 


3910 



a Present work. 



